home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / COSY_PAK / SignalProcessing / Analog / LaPlace.m < prev    next >
Encoding:
Text File  |  1992-06-08  |  42.5 KB  |  1,252 lines
  1. (*  :Title:    Forward Laplace Transform Package  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*  :Summary:    *)
  6.  
  7. (*  :Context:    SignalProcessing`Analog`LaPlace`  *)
  8.  
  9. (*  :PackageVersion:  2.6    *)
  10.  
  11. (*
  12.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  13.         Georgia Tech Research Corporation
  14.  
  15.     Permission to use, copy, modify, and distribute this software
  16.     and its documentation for any purpose and without fee is
  17.     hereby granted, provided that the above copyright notice
  18.     appear in all copies and that both that copyright notice and
  19.     this permission notice appear in supporting documentation,
  20.     and that the name of the Georgia Tech Research Corporation,
  21.     Georgia Tech, or Georgia Institute of Technology not be used
  22.     in advertising or publicity pertaining to distribution of the
  23.     software without specific, written prior permission.  Georgia
  24.     Tech makes no representations about the suitability of this
  25.     software for any purpose.  It is provided "as is" without
  26.     express or implied warranty.
  27.  *)
  28.  
  29. (*
  30.    :History:    began --  February, 1990 (adapted from "ZTransform.m")
  31.         added Dialogue ability --  February, 1991
  32. *)
  33.  
  34. (*  :Keywords:    *)
  35.  
  36. (*
  37.     :Source:    Muth, E. J.  {Transform Methods}.  1977.
  38.         Franklin, G. F. {Feedback Control of Dynamic Systems}. 1986.
  39.         Bracewell. {The Fourier Transform and Its Applications}. 1978.
  40.         Nilsson, J. W. {Electric Circuits}
  41.         Churchill, R. {Operational Mathematics}. 1958.
  42.         Oberhettinger, F., and Badii, L.  {Tables of Laplace
  43.             Transforms}, 1973.
  44.  *)
  45.  
  46. (*
  47.     :Warning:    Keep conditional clauses (code following /;) on the same line;
  48.           breaking a clause into lines separated by RETURNS will
  49.           confuse Mathematica (e.g., the active context for this
  50.           package will not be reset to Global).
  51.  *)
  52.  
  53. (*  :Mathematica Version:  1.2 or 2.0  *)
  54.  
  55. (*  :Limitation:  *)
  56.  
  57. (*
  58.     :Discussion:  All new functions have usage information.
  59.           The 1-D rule base has a total of 102 rules:
  60.  
  61.             I.   multidimensional hooks          3 rules
  62.             II.  rational transform pairs     13 rules
  63.             III. non-rational transform pairs     40 rules
  64.             IV.  transform properties         23 rules *
  65.             V.   transforms of SP structures     14 rules **
  66.             VI.  transform strategies          9 rules
  67.  
  68.           *  one rule has been commented out
  69.           ** two rules have been commented out
  70.  
  71.     The only purpose to the multidimensional hook rules is to track
  72.     the region of convergence for non-separable functions.  Without
  73.     these rules, the  z-transform  rule  base  would still properly
  74.     compute the z-transform function.   Otherwise, the CTFTransform
  75.     would  not  work  properly  for functions  like  Exp[-(6/7)^t1]
  76.     Exp[-(9/10)^t2] Delta[t1 - t2] CStep[t1, t2]  with  respect  to
  77.     variables t1 and t2.
  78.  
  79.     At each  step  in the Laplace rule base, the current expression
  80.     has a local state associated with it.  This  state  consists of
  81.     a list of six boolean values.  Each boolean value is associated
  82.     with a strategy.  If an element is True, then that strategy has
  83.     not been tried yet;  if False,  then  that strategy has already
  84.     been tried, and it will not be tried again.   Thus, local state
  85.     is used to prevent infinite loops  which would be caused by the
  86.     repetitive application  of  certain  strategy  rules.   See the
  87.     section  S T A T E    D E F I N I T I O N  below and section VI
  88.     of the rules.
  89.  
  90.     MyLaPlace[ f[t], t, s, state, tlist, slist ] is one-dimensional
  91.     Laplace transform rule base.   It returns the transform of f[t]
  92.     as a three-element list:  { F(s), Rminus, Rplus }.   The "time"
  93.     variable being  transform  is t, and s is the Laplace variable.
  94.     The  complete  collection  of  time  variables  is in tlist and
  95.     Laplace variables in slist.   Note that LaPlace drives the rule
  96.     base defined by MyLaPlace  when  computing  the forward LaPlace
  97.     transform.
  98.  *)
  99.  
  100. (*  :Functions:  LaPlace  *)
  101.  
  102.  
  103. If [ TrueQ[ $VersionNumber >= 2.0 ],
  104.      Off[ General::spell ];
  105.      Off[ General::spell1 ] ];
  106.  
  107.  
  108. (*  B E G I N     P A C K A G E  *)
  109.  
  110. BeginPackage[ "SignalProcessing`Analog`LaPlace`",
  111.           "SignalProcessing`Analog`LSupport`",
  112.           "SignalProcessing`Support`TransSupport`",
  113.           "SignalProcessing`Support`ROC`",
  114.           "SignalProcessing`Support`SigProc`",
  115.           "SignalProcessing`Support`SupCode`",
  116.           "SignalProcessing`Support`FilterSupport`" ]
  117.  
  118.  
  119. (*  U S A G E     I N F O R M A T I O N  *)
  120.  
  121. LaPlace::usage =
  122.     "LaPlace[e, t] or LaPlace[e, t, s] gives the two-sided Laplace \
  123.     transform of the expression e, which is a function of t, by \
  124.     returning an object of four slots tagged by LTransData: \
  125.     <transform>, <rminus>, <rplus>, <laplace_variables>. \
  126.     The Region of Convergence (ROC) is defined as \
  127.     <rminus> < Re{s} < <rplus>. \
  128.     Note that the returned ROC is either the actual ROC or a subset \
  129.     of the actual ROC. \
  130.     In two dimensions, LaPlace[e, {t1, t2}, {s1, s2}] \
  131.     is the same as LaPlace [ LaPlace[e, t1, s1], t2, s2 ]. \
  132.     This notation extends naturally to higher dimensions. \
  133.     Note that the right-sided transform is specified by \
  134.     multiplying the expression by CStep[t]. \
  135.     Also, LaPlaceTransform is an alias for LaPlace."
  136.  
  137. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  138.  
  139.  
  140. Begin["`Private`"]
  141.  
  142.  
  143. (*  U S E R     I N T E R F A C E  *)
  144.  
  145. (*  L operator  *)
  146. Unprotect[L]
  147. L/: TheFunction[ L[t_, s_][f_] ] := LaPlace[f, t, s]
  148. Protect[L]
  149.  
  150. (*  LaPlace  *)
  151. LaPlace/: Options[ LaPlace ] :=
  152.     { Definition -> False, Dialogue -> True,
  153.       Simplify -> True, TransformLookup -> {} }
  154.  
  155. LaPlace[e_] :=
  156.     laplacedriver[ Options[LaPlace], e ]
  157. LaPlace[e_, t_] :=
  158.     laplacedriver[ Options[LaPlace], e, t ]
  159. LaPlace[e_, t_, s_] :=
  160.     laplacedriver[ Options[LaPlace], e, t, s ]
  161. LaPlace[e_, t_, s_, op__] :=
  162.     laplacedriver[ ToList[op] ~Join~ Options[LaPlace], e, t, s ]
  163.  
  164. (*  Interface support for LaPlace  *)
  165. laplacedriver[op_, f_, args__] :=
  166.     cleanup [ laplace[op, ToContinuous[f], args],
  167.           Replace[Dialogue, op],
  168.           TrueQ[Replace[Simplify, op]] ]
  169.  
  170. laplace[op_, f_, args___] :=
  171.     Block [    {},
  172.         Replace [ ltrans[op, f, args], LaPlaceInterfaceRules ] ]
  173.  
  174. LaPlaceInterfaceRules = {
  175.     ltrans[op_, e_] :>
  176.         Message[Transform::novariables, "t (time)", GetVariables[e]],
  177.  
  178.     ltrans[op_, e_, t_] :> e    /; LTransformQ[e],
  179.     ltrans[op_, e_, t_] :>
  180.         laplace[op, e, t, DummyVariables[Length[t], Global`s] ] /;
  181.         validvarQ[t],
  182.  
  183.     ltrans[op_, e_, t_, rest___] :>
  184.         Message[ Transform::badvar, "time", LaPlace, t ] /;
  185.         ! validvarQ[t],
  186.     ltrans[op_, e_, t_, s_, rest___] :>
  187.         Message[ Transform::badvar, "s", LaPlace, s ] /;
  188.         ! validvarQ[s],
  189.  
  190.     ltrans[op_, e_, t_List, sl_] :>
  191.         MultiDTransform[ MakeLObject, laplace,
  192.                  LTransformQ, e, t, s, sl, op ],
  193.  
  194.     ltrans[op_, e_, t_Symbol, s_Symbol] :>
  195.         laplace[op, e, t, s, {t}, {s}],
  196.     ltrans[op_, e_, t_Symbol, s_Symbol, tlist_, slist_] :>
  197.         If [ SameQ[ToList[LVariables[e]], slist],
  198.              e,
  199.              LMultiDROC[laplace[op, TheFunction[e], t, s, tlist, slist],
  200.                 e] ] /;
  201.         LTransformQ[e],
  202.     ltrans[op_, e_, t_Symbol, s_Symbol, tlist_, slist_] :>
  203.         MakeLObject [ MyLaPlace [ e, t, s, initLstate[],
  204.                       tlist, slist, op ],
  205.                   s ] /;
  206.         ! LTransformQ[e]
  207. }
  208.  
  209. cleanup[ Null, dialogue_, simplify_ ] := Null
  210.  
  211. cleanup[trans_, dialogue_, simplify_ ] :=
  212.     Block [    {dialflag, retval},
  213.         retval = trans;
  214.         dialflag = SameQ[dialogue, All];
  215.         If [ simplify,
  216.              Off[Replace::rep];
  217.              retval = SPSimplify[trans, Dialogue -> dialflag];
  218.              On[Replace::rep] ];
  219.         If [ dialogue || dialflag,
  220.              Scan[explain, retval, Infinity] ];
  221.         retval ]
  222.  
  223. explain[ mylaplace[ f_, t_, rest__ ] ] :=
  224.     Message[ Transform::incomplete, "forward Laplace transform", f, t ]
  225.  
  226. validvarQ[ x_Symbol ] := True
  227. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  228. validvarQ[ x_ ] := False
  229.  
  230.  
  231. (*  S U P P O R T I N G     R O U T I N E S    F O R     R U L E     B A S E  *)
  232.  
  233. absDialogue[f_, t_, options_] :=
  234.     Block [    {dialogue, result},
  235.         dialogue = InformUserQ[ options ];
  236.         result = (f /. Abs[a_. t] :> Abs[a] t) CStep[t] +
  237.              (f /. Abs[b_. t] :> - Abs[b] t) CStep[-t];
  238.  
  239.         If [ dialogue,
  240.              Print[ "( after rewriting the two-sided expression" ];
  241.              Print[ "  ", f ];
  242.              Print[ "  as a left-sided plus a right-sided function:"];
  243.              Print[ "  ", result, " . )" ] ];
  244.  
  245.         result ]
  246.  
  247. addL[trans1_, trans2_] := AddT[LForm, trans1, trans2, -Infinity]
  248.  
  249. antiCausalL[trans_, s_] :=
  250.     Transform[ TheFunction[trans] /. s -> -s,
  251.            -GetRPlus[trans], -GetRMinus[trans] ] /;
  252.     LForm[trans]
  253.  
  254. conjL[trans_, s_] := ConjT[LForm, trans, s]
  255.  
  256. convolveL[convop_, trans1_, trans2_] :=
  257.     ConvolveT[LForm, convop, trans1, trans2]
  258.  
  259. definitionDialogue[ oldexpr_, trans_, operator_, options_ ] :=
  260.     Block [    {},
  261.         If [ InformUserQ[options],
  262.              Print[ "( after using ", operator,
  263.                 "  to apply the definition to" ];
  264.              Print[ "  ", oldexpr ];
  265.              Print[ "  to find its transform"];
  266.              Print[ "  ", trans, " )" ] ];
  267.         trans ]
  268.  
  269. integrateL[trans_, s_] := IntegrateT[LForm, trans, s, s, Infinity]
  270.  
  271. lDerivative[trans_, s_, m_] := scaleL[DerivativeT[LForm, trans, s, m], (-1)^m]
  272.  
  273. lineImpulseMDL[trans_, s_, slist_, sleft_] :=
  274.     LineImpulsemDT[LForm, trans, s, slist, sleft, Plus]
  275.  
  276. lMultiplyByExp[trans_, s_, a_] :=
  277.     Block [    {newrm, newrp, transform},
  278.         newrm = GetRMinus[trans] - Re[a];
  279.         newrp = GetRPlus[trans] - Re[a];
  280.         transform = Transform[ TheFunction[trans], newrm, newrp ];
  281.         substituteForL[transform, s, s+a] ] /;
  282.     LForm[trans]
  283.  
  284. multL[trans1_, trans2_] := MultT[LForm, trans1, trans2, -Infinity]
  285.  
  286. scaleL[trans_, c_] := ScaleT[LForm, trans, c]
  287.  
  288. similarityL[x_, s_, scale_] :=
  289.     Block [    {adjrm, adjrp, newx},
  290.         adjrm = GetRMinus[x];
  291.         If [ ! InfinityQ[adjrm], adjrm *= scale ];
  292.         adjrp = GetRPlus[x];
  293.         If [ ! InfinityQ[adjrp], adjrp *= scale ];
  294.         newx = Transform[ TheFunction[x],
  295.                   Min[adjrm, adjrp],
  296.                   Max[adjrm, adjrp] ];
  297.         scaleL[ substituteForL[ newx, s, s/scale ], 1/Abs[scale] ] ] /;
  298.     LForm[x]
  299.  
  300. similarityQ[f_, t_] :=
  301.     Block [    {scale},
  302.         scale = ScalingFactor[f, t];
  303.         (! SameQ[scale, 1]) && (! SameQ[scale, 0]) ]
  304.  
  305. subL[trans1_, trans2_] := SubT[LForm, trans1, trans2, -Infinity]
  306.  
  307. substituteForL[trans_, s_, news_] := SubstituteForT[LForm, trans, s, news]
  308.  
  309. summationL[trans_, op_ ] :=
  310.     Block [    {fun, rm, rp},
  311.         fun = TheFunction[trans];
  312.         rm = GetRMinus[trans];
  313.         rp = GetRPlus[trans];
  314.         Transform[ op[fun], rm, rp ] ] /;
  315.     LForm[trans]
  316.  
  317. (* one-sided function *)
  318. timeDerivative[x_, s_, m_, f_, t_] :=
  319.     Block [    {curderivative, i, result, rm, rp, xfun},
  320.         rm = Max[0, GetRMinus[x]];
  321.         rp = GetRPlus[x];
  322.         xfun = TheFunction[x];
  323.  
  324.         result = s^m xfun  -  s^(m-1) Limit[f, t -> 0];
  325.         curderivative = f;
  326.         For [ i = 2, i <= m, i++,
  327.               curderivative = D[curderivative, t];
  328.               result -= s^(m-i) Limit[curderivative, t -> 0] ];
  329.         Transform[result, rm, rp] ] /;
  330.     LForm[x] && IntegerQ[m]
  331.  
  332.  
  333. (*  S T A T E     D E F I N T I O N  *)
  334.  
  335. (*  In order of appearance in rule base  *)
  336.  
  337. polyfactorfield = 1        (* factor denominator if denom. is rat. poly  *)
  338. partialfractionsfield = 2    (* partial fractions if denom. is rat. poly.  *)
  339. expandfield = 3            (* apply Expand[] to expression              *)
  340. expandallfield = 4        (* apply ExpandAll[] to expression          *)
  341. collectallfield = 5        (* apply Collect to every level of expression *)
  342. stepfield = 6            (* multiply by (CStep[t] + CStep[-t])          *)
  343. thefunfield = 7            (* apply TheFunction[] to expression          *)
  344. definitionfield = 8        (* apply definition of Laplace transform      *)
  345.  
  346. statevariables = 8
  347.  
  348. initLstate[] := Table[True, {statevariables}]
  349. nullLstate[] := Table[False, {statevariables}]
  350.  
  351.  
  352. (*   Driver for one-dimensional rule base            *)
  353. (*   First, convert all z /(z - a) forms to 1/(1 - a z^-1)    *) 
  354. (*    since the rule base favors terms with z^-1 terms    *)
  355. (*   Loop until the expression to be inverted does not change.    *)
  356. (*   Use fixUp to remove the TransformLookup option before    *)
  357. (*    the options are passed on to the rule base; otherwise,    *)
  358. (*    the rule base could loop indefinitely.            *)
  359. MyLaPlace[ f_, t_, s_, st_, tlist_, slist_, op_ ] :=
  360.     Block [ {laste = Null, newe, newf, newlrules, trace},
  361.  
  362.         (* generate the forward Laplace transform rules *)
  363.         newlrules = TransformFixUp[ t, s, op, mylaplace, True,
  364.                         LaPlace, -Infinity, Infinity ];
  365.         LaPlaceRules = Join[ newlrules, OriginalLaPlaceRules ];
  366.  
  367.         (* determine the level of dialogue *)
  368.         trace = SameQ[ Replace[Dialogue, op], All ];
  369.  
  370.         (* separate Sqrt[b t] forms into Sqrt[b] Sqrt[t] forms *)
  371.         newf = f /. ( Sqrt[b_ t] :> Sqrt[b] Sqrt[t] /; FreeQ[b, t] );
  372.  
  373.         (* take the 1-D Laplace transform *)
  374.         newe = mylaplace[newf, t, s, st, tlist, slist, fixUp[op]];
  375.         While [ ! SameQ[ laste, newe ],
  376.             If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
  377.             laste = newe;
  378.             newe = MapAll[transform, laste] ];
  379.         If [ trace, Print[ newe ] ];
  380.  
  381.         newe ]
  382.  
  383. (*  fixUp --  remove redundant options and the option TransformLookup  *)
  384. fixUp[ op_ ] := { Definition -> Replace[Definition, op],
  385.           Dialogue -> Replace[Dialogue, op],
  386.           Simplify -> Replace[Simplify, op] }
  387.  
  388. transform[ expr_ ] :=
  389.     If [ SameQ[Head[expr], mylaplace],
  390.          Replace [ expr, LaPlaceRules ],
  391.          expr ]
  392.  
  393.  
  394. (* Format intermediate forms so that output from Dialogue -> All is readable *)
  395.  
  396. Format[ mylaplace[ x_, t_, s_, st_, tlist_, slist_, op_ ] ] := 
  397.     SequenceForm[ ColumnForm[{"L",
  398.                   "  " ~StringJoin~ ToString[t]}],
  399.               { x } ]
  400.  
  401. Format[ antiCausalL[trans_, s_] ] :=
  402.     SequenceForm[ {trans}, Subscript[s -> -s], Subscript[" and flip ROC "] ]
  403.  
  404. Format[ lMultiplyByExp[trans_, s_, a_] ] :=
  405.     SequenceForm[ {trans},
  406.               Subscript[s -> s + a],
  407.               Subscript[" and shift ROC "] ]
  408.  
  409. Format[ similarityL[trans_, s_, scale_] ] :=
  410.     SequenceForm[ {trans},
  411.               Subscript[s -> s/scale],
  412.               Subscript[" and scale transform & ROC "] ]
  413.  
  414. tdfmt = "`` - Summation[i, 0, ``, 1][`` Limit[ D[``, {``, i}], `` -> 0 ] ]"
  415.  
  416. Format[ timeDerivative[x_, s_, m_, f_, t_] ] :=
  417.     ToString[ StringForm[ tdfmt, s^m x, m - 1, s^(m - "i"), f, t, t ] ]
  418.  
  419.  
  420. (*  B E G I N     R U L E     B A S E  *)
  421.  
  422.  
  423. LaPlaceRules = {}
  424.  
  425.  
  426. OriginalLaPlaceRules = {
  427.  
  428.  
  429. (*  I.    M U L T I D I M E N S I O N A L     H O O K S  *)
  430.  
  431.  
  432. (*        A.  Region of convergence specification for dimension given    *)
  433. (*        by the variable n.  Allows multi-dimensional LaPlace    *)
  434. (*        transforms like a^t1 Delta[t1 - t2] CStep[t1,t2] to be    *)
  435. (*        taken --  see Line impulses                *)
  436. mylaplace[ SignalProcessing`ROCinfo[t_, rm_:-Infinity , rp_:Infinity], t_, s_,
  437.        st_, tlist_, slist_, options_ ] :>
  438.     Transform[ 1, rm, rp ],
  439.  
  440.  
  441. (*      B.  Multidimensional line Deltas become LineImpulses.        *)
  442. mylaplace[ f_. Delta[k_. t1_Symbol + l_. t2_Symbol], t_, s_,
  443.        st_, tlist_, slist_, options_ ] :>
  444.     mylaplace[ f LineImpulse[{t1,t2}, {k,-l}], t, s,
  445.            st, tlist, slist, options ] /;
  446.     SubsetQ[{t1, t2}, tlist] && MemberQ[{t1, t2}, t],
  447.  
  448.  
  449. (*      C.  Line impulses.                        *)
  450. mylaplace[ f_. LineImpulse[varlist_, coefflist_], t_, s_,
  451.        st_, tlist_, slist_, options_ ] :>
  452.     Block [    {functionoft, slistmD},
  453.         functionoft = f /. ReplaceWith[varlist, t / coefflist];
  454.         slistmD = AssociateItem[varlist, tlist, slist];
  455.         lineImpulseMDL [ mylaplace [ functionoft, t, s,
  456.                          st, tlist, slist, options ],
  457.                      s, slistmD, Complement[varlist,{t}]] ] /;
  458.     FreeQ[f, LineImpulse[a__]],
  459.  
  460.  
  461.  
  462.  
  463. (*  II.   R A T I O N A L     L A P L A C E     T R A N S F O R M S    *)
  464.  
  465.  
  466. (*      A.  Lookup rules for the zero, step, and impulse functions.    *)
  467. mylaplace[ 0, t_, s_, st_, tlist_, slist_, options_ ] :>
  468.     Transform[ 0, -Infinity, Infinity ],
  469.  
  470. mylaplace[ CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  471.     Transform[ 1/s, 0, Infinity ],
  472.  
  473. mylaplace[ f_. Delta[a_. t_ + t0_.], t_, s_, st_, tlist_, slist_, options_ ] :>
  474.     Transform[ (f /. t -> -t0/a) Exp[s t0/a] / Abs[a],
  475.             -Infinity, Infinity ] /;
  476.     FreeQ[{a,t0}, t],
  477.  
  478. mylaplace[ Unit[n_][t_ + t0_.], t_, s_, st_, tlist_, slist_, options_ ] :>
  479.     Transform[ s^n Exp[s t0], -Infinity, Infinity ] /;
  480.     IntegerQ[n] && n > 0 && FreeQ[t0, t],
  481.  
  482. mylaplace[ t_^n_. CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  483.     Block [    {dialogue},
  484.         dialogue = SameQ[Replace[Dialogue, options], All];
  485.         Assuming[ n > 0, dialogue ];
  486.         If [ dialogue && ! IntegerQ[n],
  487.              Print[ "and that ", n, " is an integer" ] ];
  488.         Transform[ n! / s^(n+1), 0, Infinity ] ] /;
  489.     Implies[ NumberQ[n], IntegerQ[n] && n > 0 ],
  490.  
  491.  
  492. (*      B.  Sinusoidal functions --  note that cross terms like    *)
  493. (*        cos(at) cos(bt) are handled by properties        *)
  494. mylaplace[ Sin[b_. + w_. t_] CStep[t_], t_, s_,
  495.        st_, tlist_, slist_, options_ ] :>
  496.     Transform[ (s Sin[b] + w Cos[b]) / (s^2 + w^2), 0, Infinity ] /;
  497.      FreeQ[{b,w}, t],             (*  Modified from [Franklin, 591]  *)
  498.  
  499. mylaplace[ Cos[b_. + w_. t_] CStep[t_], t_, s_,
  500.        st_, tlist_, slist_, options_ ] :>
  501.     Transform[ (s Cos[b] - w Sin[b]) / (s^2 + w^2), 0, Infinity ] /;
  502.      FreeQ[{b,w}, t],             (*  Modified from [Franklin, 591]  *)
  503.  
  504. mylaplace[ a_. Sinh[b_. + w_. t_] CStep[t_], t_, s_,
  505.        st_, tlist_, slist_, options_ ] :>
  506.     Transform[a (s Sinh[b] + w Cosh[b]) / (s^2 - w^2), Abs[w], Infinity] /;
  507.      FreeQ[{a,b,w}, t],             (*  Modified from [Franklin, 591]  *)
  508.  
  509. mylaplace[ a_. Cosh[b_. + w_. t_] CStep[t_], t_, s_,
  510.        st_, tlist_, slist_, options_ ] :>
  511.     Transform[a (s Cosh[b] - w Sinh[b]) / (s^2 - w^2), Abs[w], Infinity] /;
  512.      FreeQ[{a,b,w}, t],             (*  Modified from [Franklin, 591]  *)
  513.  
  514. mylaplace[ Sin[a_. t_]^2 CStep[t_], t_, s_,
  515.        st_, tlist_, slist_, options_ ] :>
  516.     Transform[ 2 a^2 / ( s (s^2 + 4 a^2) ), 2 Abs[Im[a]], Infinity ] /;
  517.     FreeQ[a, t],                 (*  [O & B, Formula 7.38]  *)
  518.  
  519. mylaplace[ Cos[a_. t_]^2 CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  520.     Transform[ (s^2 + 2 a^2)/(s (s^2 + 4 a^2)), 2 Abs[Im[a]], Infinity ] /;
  521.     FreeQ[a, t],                 (*  [O & B, Formula 7.39]  *)
  522.  
  523. (*      C.  Constant function                        *)
  524. mylaplace[ a_, t_, s_, st_, tlist_, slist_, options_ ] :>
  525.     Transform[ a Delta[s], 0, 0 ] /;
  526.     FreeQ[a, t],
  527.  
  528. (*      D.  LaGuerre polynomials                    *)
  529. mylaplace[ LaguerreL[n_, t_ + b_.] CStep[t_], t_, s_,
  530.         st_, tlist_, slist_, options_ ] :>
  531.     Transform[ Exp[b s] ((s - 1)/s)^n / s, 0, Infinity ] /;
  532.     FreeQ[{b,n}, t],
  533.  
  534. (*      E.  Other rational transform pairs                *)
  535. mylaplace[ t_^j_. BesselJ[j_, a_. t_] CStep[t_], t_, s_,
  536.         st_, tlist_, slist_, options_ ] :>
  537.     Transform[ (2 a)^j Gamma[j + 1/2] / ( Sqrt[Pi] (s^2 + a^2)^(j + 1/2) ),
  538.            0, Infinity ] /;
  539.     FreeQ[{a,j}, t],
  540.  
  541.  
  542.  
  543.  
  544. (*  III.  N O N - R A T I O N A L     L A P L A C E     T R A N S F O R M S  *)
  545.  
  546.  
  547. (*      A.  "Time" functions whose denominators are rational polys    *)
  548.  
  549. (*          1.  Put functions with rational polynomial denominators    *)
  550. (*          in standard form (factor denominator, expand numer.)    *)
  551. mylaplace[ p_, t_, s_, st_, tlist_, slist_, options_ ] :>
  552.     mylaplace [ Expand[Numerator[p]] / ( Factor[Denominator[p]] ), t, s,
  553.             SetStateField[st, polyfactorfield, False],
  554.             tlist, slist, options ] /;
  555.     GetStateField[st, polyfactorfield] &&
  556.       PolynomialQ[Denominator[p], t] && RationalFunctionQ[p, t],
  557.  
  558. (*          2.  Normalize the factors of the denominator        *)
  559. (*          a.  for numeric k, numeric j                *)
  560. mylaplace[ ((a_ + b_ t_^j_.)^k_) x_., t_, s_,
  561.         st_, tlist_, slist_, options_ ] :>
  562.     b^k mylaplace[ x ( a/b + t^j )^k, t, s, st, tlist, slist, options ] /;
  563.     FreeQ[{a,b},t] && IntegerQ[k] && (k < 0) && IntegerQ[j] && (j > 0),
  564.  
  565. (*          b.  for symbolic k, numeric j                *)
  566. mylaplace[ x_. / ((a_ + b_ t_^j_.)^k_.), t_, s_,
  567.        st_, tlist_, slist_, options_ ] :>
  568.     b^(- k) mylaplace[x / (a/b + t^j)^k, t, s, st, tlist, slist, options] /;
  569.     FreeQ[{a,b,k},t] && IntegerQ[j] && (j > 0),
  570.  
  571. (*          3.  Partial fractions decomposition            *)
  572. mylaplace[ p_, t_, s_, st_, tlist_, slist_, options_ ] :>
  573.     mylaplace [ Apart[p, t], t, s, 
  574.             SetStateField[st, partialfractionsfield, False],
  575.             tlist, slist, options ] /;
  576.     GetStateField[st, partialfractionsfield] &&
  577.       PolynomialQ[Denominator[p], t] && RationalFunctionQ[p, t],
  578.  
  579. (*          4.  Denominator sections                    *)
  580. mylaplace[ CStep[t_] / ( t_ + a_ ), t_, s_, st_, tlist_, slist_, options_ ] :>
  581.     Transform[ - Exp[a s] ExpIntegralEi[- a s], 0, Infinity ] /;
  582.     FreeQ[a, t],                  (*  [O & B, Formula 2.9]  *)
  583.  
  584. mylaplace[ CStep[t_] / ( t_ + a_ )^2, t_, s_,
  585.         st_, tlist_, slist_, options_ ] :>
  586.     Transform[ 1 / a + s Exp[a s] ExpIntegralEi[- a s], 0, Infinity ] /;
  587.     FreeQ[a, t],
  588.  
  589. mylaplace[ CStep[t_] / ( t_^2 + b_ ), t_, s_,
  590.         st_, tlist_, slist_, options_ ] :>
  591.     Block [    {a, trans},
  592.         a = Sqrt[b];
  593.         trans = CosIntegral[a s] Sin[a s] / a -
  594.             SinIntegral[a s] Cos[a s] / a;
  595.         Transform[trans, 0, Infinity] ] /;
  596.     FreeQ[b, t] && Positive[b],         (*  [O & B, Formula 2.16]  *)
  597.  
  598. mylaplace[ CStep[t_] / ( t_^2 + b_ ), t_, s_,
  599.         st_, tlist_, slist_, options_ ] :>
  600.     Block [    {g, trans},
  601.         a = Sqrt[- b];
  602.         trans = ExpIntegralEi[- a s] Exp[a s] / ( 2 a ) -
  603.             (I Pi + ExpIntegralEi[a s] Exp[- a s] ) / ( 2 a );
  604.         Transform[trans, 0, Infinity] ] /;
  605.     FreeQ[{a,b}, t] && Negative[b],         (*  [O & B, Formula 2.21]  *)
  606.  
  607. mylaplace[ (a_. t_ + b_.) CStep[t_] / (t_^2 + d_), t_, s_,
  608.        st_, tlist_, slist_, options_ ] :>
  609.     Block [ {c, trans},
  610.         c = Sqrt[d];
  611.         trans = - CosIntegral[c s] (a Cos[c s] - b Sin[c s] / c) -
  612.             SinIntegral[c s] (a Sin[c s] + b Cos[c s] / c);
  613.         Transform[trans, 0, Infinity] ] /;
  614.     FreeQ[{a,b,c,d}, t] && Positive[d],  (*  [O & B, Formula 2.19]  *)
  615.  
  616. mylaplace[ (a_. t_ + b_.) CStep[t_] / (t_^2 + d_), t_, s_,
  617.        st_, tlist_, slist_, options_ ] :>
  618.     Block [    {c, trans},
  619.         c = Sqrt[- d];
  620.         trans = - ExpIntegralEi[- c s] Exp[c s] (a - b/c) / (2 f) -
  621.             ExpIntegralEi[c s] Exp[- c s] (a + b/c) / (2 f );
  622.         Transform[trans, 0, Infinity] ] /;
  623.     FreeQ[{a,b,d}, t] && Negative[d],    (*  [O & B, Formula 2.23]  *)
  624.  
  625. mylaplace[ CStep[t_] / t_, t_, s_, st_, tlist_, slist_, options_ ] :>
  626.     Transform[ - Log[s], 0, Infinity ],      (*  [Churchill, 324]  *)
  627.  
  628. mylaplace[ (t_ + b_.)^k_. CStep[t_ + b_.], t_, s_,
  629.         st_, tlist_, slist_, options_ ] :>
  630.     Block [    {},
  631.         Assuming[ Positive[k + 1], options];
  632.         Transform[ Exp[b s] Gamma[k+1] / s^(k+1), 0, Infinity ] ] /;
  633.     FreeQ[b, t] && Implies[NumberQ[N[k]], N[k > -1]],
  634.                           (*  [O & B, Formula 3.3]  *)
  635.  
  636.  
  637. (*      B.  "Time" functions whose denominators contains powers of t    *)
  638.  
  639. mylaplace[ CPulse[l_, t_ + k_.] / Sqrt[ t_ ( t_ + l_ ) ],
  640.        t_, s_, st_, tlist_, slist_, options_ ] :>
  641.     Transform[ -I Pi Exp[s k] Exp[l s / 2] BesselI[0, l s / 2],
  642.             -Infinity, Infinity ] /;
  643.     FreeQ[{k,l}, t],            (*  [Churchill, 323 & 331]  *)
  644.  
  645. mylaplace[ t_ CPulse[l_, t_ + k_.] / Sqrt[ t_ ( t_ + l_ ) ], t_, s_,
  646.        st_, tlist_, slist_, options_ ] :>
  647.     Transform[ I Pi l Exp[s k] Exp[l s / 2] BesselI[0, l s / 2] -
  648.             I Pi l Exp[l s / 2] BesselI[1, l s / 2],
  649.             -Infinity, Infinity ] /;
  650.     FreeQ[{k,l}, t],            (*  [Churchill, 323 & 331]  *)
  651.  
  652.  
  653. (*      C.  Gaussian function                        *)
  654. mylaplace[ Exp[a_. t_^2] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  655.     Transform[ Sqrt[Pi / -a] Exp[- s^2 / (4 a)] *
  656.             ( 1 - Erf[s / (2 Sqrt[-a])] ), -Infinity, Infinity ] /;
  657.     FreeQ[a, t],                  (*  [O & B, Formula 5.41] *)
  658.  
  659.  
  660. (*      D.  Common waves                        *)
  661.  
  662. (*          1.  square waves                        *)
  663. mylaplace[ SquareWave[c_, t_] CStep[t_], t_, s_,
  664.        st_, tlist_, slist_, options_ ] :>
  665.     Transform[Tanh[c s / 2] / s, 0, Infinity] /;
  666.     FreeQ[c, t],                   (*  [Churchill, 54]  *)
  667.  
  668. (*          2.  full-wave rectified sinusoids                *)
  669. mylaplace[ Abs[Sin[a_. t_]] CStep[t_], t_, s_,
  670.        st_, tlist_, slist_, options_ ] :>
  671.     Transform[a Coth[Pi s / (2 a)] / (s^2 + a^2), Abs[Im[a]], Infinity] /;
  672.     FreeQ[a, t],                  (*  [O & B, Formula 7.3]  *)
  673.  
  674. mylaplace[ Abs[Cos[a_. t_]] CStep[t_], t_, s_,
  675.        st_, tlist_, slist_, options_ ] :>
  676.     Transform[ ( s + a Csch[Pi s / (2 a)] ) / (s^2 + a^2),
  677.             Abs[Im[a]], Infinity ] /;
  678.     FreeQ[a, t],                  (*  [O & B, Formula 7.4]  *)
  679.  
  680.  
  681. (*      E.  Bessel functions                        *)
  682.  
  683. (*  Bessel functions [Churchill, 327] and [O & B, Formulas 14.1/15.1]    *)
  684. (*  Note that formula 14.1 in O & B has a typo in it: -v should be v    *)
  685.  
  686. mylaplace[ BesselI[v_, a_. t_] CStep[t_], t_, s_,
  687.        st_, tlist_, slist_, options_ ] :>
  688.     Block [    {},
  689.         Assuming[ Positive[Re[v] + 1], options ];
  690.         Transform[ (s - Sqrt[s^2 - a^2])^v / (a^v Sqrt[s^2 - a^2]),
  691.                Abs[Re[a]], Infinity ] ] /;
  692.     FreeQ[{a,v}, t] && Implies[NumberQ[N[v]], N[Re[v] > -1]],
  693.  
  694. mylaplace[ BesselJ[v_, a_. t_] CStep[t_], t_, s_,
  695.        st_, tlist_, slist_, options_ ] :>
  696.     Block [    {},
  697.         Assuming[ Positive[Re[v] + 1], options];
  698.         Transform[ (Sqrt[s^2 + a^2] - s)^v / (a^v Sqrt[s^2 + a^2]),
  699.                 Abs[Im[a]], Infinity ] ] /;
  700.     FreeQ[{a,v}, t] && Implies[NumberQ[N[v]], N[Re[v] > -1]],
  701.  
  702. mylaplace[ BesselJ[v_, a_. t_] CStep[t_] / t_, t_, s_,
  703.        st_, tlist_, slist_, options_ ] :>
  704.     Block [    {},
  705.         Assuming[ Positive[Re[v]], options ];
  706.         Transform[ ( Sqrt[s^2 + a^2] - s )^v / ( v a^v ),
  707.                 Abs[Im[a]], Infinity ] ] /;
  708.     FreeQ[{a,v}, t] && Implies[NumberQ[N[v]], Re[N[v]] > 0],
  709.  
  710. mylaplace[ BesselK[0, j_. Sqrt[t_]] CStep[t_] / Sqrt[t_], t_, s_,
  711.        st_, tlist_, slist_, options_ ] :>
  712.     Transform[ Sqrt[Pi] Exp[j^2 / (8 s)] BesselK[0, j^2 / (8 s)] /
  713.             (2 Sqrt[s]), 0, Infinity ] /;
  714.     FreeQ[j, t],
  715.  
  716.  
  717. (*      F.  Exponential forms                        *)
  718.  
  719. (*          1.  exponential integral  [Churchill, 329]        *)
  720. mylaplace[ ExpIntegralEi[a_. t_] CStep[t_], t_, s_,
  721.        st_, tlist_, slist_, options_ ] :>
  722.     Block [    {},
  723.         Assuming[ Negative[a], options ];
  724.         Transform[ - Log[ 1 - s / a ] / s, 0, Infinity ] ] /;
  725.     FreeQ[a, t] && Implies[NumberQ[N[a]], N[a < 0]],
  726.                            (* [Churchill, 329]  *)
  727.  
  728. (*          2.  error functions, erf                    *)
  729. mylaplace[ Erf[a_. Sqrt[t_]] CStep[t_], t_, s_,
  730.        st_, tlist_, slist_, options_ ] :>
  731.     Transform[ a / ( s Sqrt[s + a^2] ), 0, Infinity ] /;
  732.     FreeQ[a, t],             (*  Modified [Churchill, 326]  *)
  733.  
  734. mylaplace[ Erf[a_. t_ + b_.] CStep[t_], t_, s_,
  735.        st_, tlist_, slist_, options_ ] :>
  736.     Transform[ ( Exp[s^2 / (4 a^2) + s b / a] ( 1 - Erf[b + s / (2 a)] ) +
  737.               Erf[b] ) / s, 0, Infinity ] /;
  738.     FreeQ[{a,b}, t],            (*  [O & B, Formula 17.59]  *)
  739.  
  740. mylaplace[ Exp[a_/t_] CStep[t_], t_, s_, 
  741.        st_, tlist_, slist_, options_ ] :>
  742.     Transform[ -2 I Sqrt[a] BesselK[-1, -2 I Sqrt[a] Sqrt[s]] / Sqrt[s],
  743.            0, Infinity ] /;
  744.     FreeQ[a, t],        (* from Mathematica's Integrate command *)
  745.  
  746. mylaplace[ Exp[a_/t_] CStep[t_] / Sqrt[t_], t_, s_, 
  747.        st_, tlist_, slist_, options_ ] :>
  748.     Transform[ Exp[2 I Sqrt[a] Sqrt[s]] Sqrt[Pi] / Sqrt[s], 0, Infinity ] /;
  749.     FreeQ[a, t],        (* from Mathematica's Integrate command *)
  750.  
  751.  
  752. (*      G.  Trigonometric functions                    *)
  753.  
  754. (*          1.  sinc-like forms                    *)
  755. mylaplace[ Sin[a_. t_] CStep[t_] / t_, t_, s_,
  756.         st_, tlist_, slist_, options_ ] :>
  757.     Transform[ ArcTan[ a / s ], Abs[Im[a]], Infinity ] /;
  758.     FreeQ[a, t],                  (*  [O & B, Formula 7.5]  *)
  759.  
  760. mylaplace[ Sin[k_. Sqrt[t_]] CStep[t_] / t_, t_, s_,
  761.         st_, tlist_, slist_, options_ ] :>
  762.     Transform[ Pi Erf[ k / ( 2 Sqrt[s] ) ], -Infinity, Infinity ] /;
  763.     FreeQ[k, t],                  (*  [Churchill, 330]  *)
  764.  
  765. (*          2.  sine integral Si                    *)
  766. mylaplace[ SinIntegral[a_. t_ + b_.] CStep[t_], t_, s_,
  767.        st_, tlist_, slist_, options_ ] :>
  768.     Transform[ Exp[b s / a] ArcTan[a / s] / s, 0, Infinity ] /;
  769.     FreeQ[{a,b}, t],              (*  [Churchill, 330]  *)
  770.  
  771. (*          3.  cosine integral Ci                    *)
  772. mylaplace[ CosIntegral[a_. t_] CStep[t_], t_, s_,
  773.        st_, tlist_, slist_, options_ ] :>
  774.     Transform[ ( 2 Log[a] - Log[s^2 + a^2] ) / ( 2 s ), 0, Infinity ] /;
  775.     Positive[a],                  (*  [Churchill, 330]  *)
  776.  
  777. mylaplace[ CosIntegral[a_. t_] CStep[t_], t_, s_,
  778.        st_, tlist_, slist_, options_ ] :>
  779.     Transform[ - Log[1 + s^2/a^2] / (2 s), 0, Infinity ] /;
  780.     FreeQ[a, t],                  (*  [Churchill, 330]  *)
  781.  
  782. (*          4.  inverse trigonometric functions            *)
  783. mylaplace[ ArcTan[a_. t_] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  784.     Transform[ (CosIntegral[s/a] Sin[s/a] - SinIntegral[s/a] Cos[s/a]) / s,
  785.             0, Infinity ] /;
  786.     FreeQ[a, t],                  (*  [O & B, Formula 8.4]  *)
  787.  
  788. mylaplace[ ArcCot[a_. t_] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  789.     Transform[ ( Pi/2 - CosIntegral[s/a] Sin[s/a] -
  790.               SinIntegral[s/a] Cos[s/a] ) / s, 0, Infinity ] /;
  791.     FreeQ[a, t],                  (*  [O & B, Formula 8.5]  *)
  792.  
  793.  
  794. (*      H.  Logarithmic forms                        *)
  795.  
  796. mylaplace[ Log[a_. t_] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  797.     Transform[ - a ( EulerGamma + Log[s / a] ) / s, 0, Infinity ] /;
  798.     FreeQ[a, t],             (*  Modified [O & B, Formula 6.1]  *)
  799.  
  800. mylaplace[ Log[a_. t_]^2 CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  801.     Transform[ a ( Pi^2 / 6 + ( EulerGamma + Log[s / a] )^2 ) / s,
  802.             0, Infinity ] /;
  803.     FreeQ[a, t],                  (*  [O & B, Formula 6.2]  *)
  804.  
  805. mylaplace[ Log[a_ + b_. t_] CStep[t_], t_, s_,
  806.        st_, tlist_, slist_, options_ ] :>
  807.     Block [    {},
  808.         Assuming[ Negative[ Abs[Arg[b/a]] - Pi ], options ];
  809.         Transform[ (Log[a] - Exp[a s / b] ExpIntegralEi[- a s / b]) / s,
  810.                 0, Infinity ] ] /;
  811.     FreeQ[{a,b},t] && Implies[NumberQ[N[a]] && NumberQ[N[b]], N[Abs[Arg[b/a]] < Pi]],                      (*  [O & B, Formula 6.7]  *)
  812.  
  813. mylaplace[ Log[t_^2 + b_] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  814.     Block [    {a, trans},
  815.         Assuming[ Positive[b], options];
  816.         a = Sqrt[b];
  817.         trans = 2 ( Log[a] - CosIntegral[a s] Cos[a s] -
  818.                 SinIntegral[a s] Sin[a s] ) / s;
  819.         Transform[trans, 0, Infinity] ] /;
  820.     FreeQ[b, t] && Implies[ NumberQ[N[b]], Positive[b] ],
  821.                          (*  [O & B, Formula 6.11]  *)
  822.  
  823. mylaplace[ Log[1 + b_. t_^2] CStep[t_] / t_, t_, s_,
  824.        st_, tlist_, slist_, options_ ] :>
  825.     Block [    {a, trans},
  826.         a = Sqrt[b];
  827.         trans = CosIntegral[s/a]^2 + SinIntegral[s/a]^2;
  828.         Transform[trans, 0, Infinity] ] /;
  829.     FreeQ[b, t],                 (*  [O & B, Formula 6.15]  *)
  830.  
  831.  
  832.  
  833.  
  834.  
  835. (*  IV.   T R A N S F O R M     P R O P E R T I E S            *)
  836.  
  837.  
  838. (*      A.  pick off constants                    *)
  839. mylaplace[ c_ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  840.     scaleL[mylaplace[x, t, s, st, tlist, slist, options], c] /;
  841.     FreeQ[c, t] && ! FreeQ[x, t],
  842.  
  843.  
  844. (*      B.  additivity -- resulting ROC is intersection of two ROC's    *)
  845. mylaplace[ (x_ + y_) CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  846.     addL[mylaplace[x CStep[t], t, s, st, tlist, slist, options],
  847.          mylaplace[y CStep[t], t, s, st, tlist, slist, options]],
  848.  
  849. mylaplace[ x_ + y_, t_, s_, st_, tlist_, slist_, options_ ] :>
  850.     addL[mylaplace[x, t, s, st, tlist, slist, options],
  851.          mylaplace[y, t, s, st, tlist, slist, options]],
  852.  
  853. mylaplace[ (x_ + y_) CStep[t_] / c_, t_, s_, st_, tlist_, slist_, options_ ] :>
  854.     addL[mylaplace[x CStep[t] / c, t, s, st, tlist, slist, options],
  855.          mylaplace[y CStep[t] / c, t, s, st, tlist, slist, options]],
  856.  
  857. mylaplace[ (x_ + y_) / c_, t_, s_, st_, tlist_, slist_, options_ ] :>
  858.     addL[mylaplace[x / c, t, s, st, tlist, slist, options],
  859.          mylaplace[y / c, t, s, st, tlist, slist, options]],
  860.  
  861.  
  862. (*      C.  Time reversal                        *)
  863. mylaplace[ x_. CStep[-t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  864.     antiCausalL [ mylaplace [ (x /. t -> - t) CStep[t], t, s,
  865.                   st, tlist, slist, options ],
  866.               s ],
  867.  
  868.  
  869. (*      D.  Shift and delays --  every continuous function should    *)
  870. (*          be a function times either a CStep  or an Delta, and    *)
  871. (*          the rule base already handles the Delta case.        *)
  872. mylaplace[ f_. CStep[a_ t_ + b_.], t_, s_, st_, tlist_, slist_, options_ ] :>
  873.     mylaplace[ f CStep[t + b/a] CStep[a] + f CStep[-t - b/a] CStep[-a],
  874.            t, s, st, tlist, slist, options ] /;
  875.     FreeQ[{a, b}, t],                 (* [Bracewell, 67] *)
  876.  
  877. mylaplace[ f_. CStep[t_ + m_], t_, s_, st_, tlist_, slist_, options_ ] :>
  878.     scaleL [ mylaplace [ (f /. t -> t - m) CStep[t], t, s,
  879.                  st, tlist, slist, options ],
  880.              Exp[m s] ] /;
  881.     FreeQ[m, t],
  882.  
  883. mylaplace[ f_. CStep[m_ - t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  884.     scaleL [ mylaplace [ (f /. t -> t + m) CStep[-t], t, s,
  885.                  st, tlist, slist, options ],
  886.              Exp[- m s] ] /;
  887.     FreeQ[m, t],
  888.  
  889. mylaplace[ f_[t_ + b_], t_, s_, st_, tlist_, slist_, options_ ] :>
  890.     scaleL[ mylaplace[ f[t], t, s, st, tlist, slist, options],
  891.         Exp[b s] ] /;
  892.     FreeQ[b, t],
  893.  
  894. (*
  895. mylaplace[ f_[t_ + b_] CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  896.     scaleL[ mylaplace[ f[t] CStep[t], t, s, st, tlist, slist, options],
  897.         Exp[b s] ] /;
  898.     FreeQ[b, t],
  899.  *)
  900.  
  901.  
  902. (*      E.  Multiplication by c^(at+b) --  covers exponential case,    *)
  903. (*            b can be a function of t                    *)
  904. mylaplace[ c_^(b_. + a_. t_) x_., t_, s_, st_, tlist_, slist_, options_ ] :>
  905.     lMultiplyByExp [ mylaplace[c^b x, t, s, st, tlist, slist, options],
  906.              s, -a Log[c] ]  /;
  907.     FreeQ[{a,c}, t],           (* Modified [Franklin, 590]  *)
  908.  
  909.  
  910. (*      F.  Multiplication by t^m [Franklin, 590]            *)
  911. mylaplace[ t_^m_. x_., t_, s_, st_, tlist_, slist_, options_ ] :>
  912.     lDerivative[ mylaplace[x, t, s, st, tlist, slist, options], s, m ] /;
  913.     IntegerQ[m] && m > 0,
  914.  
  915.  
  916. (*      G.  Modulation                        *)
  917.  
  918. (*          1.  multiplication by a cosine function [Bracewell, 224]    *)
  919. mylaplace[ Cos[b_ + w_. t_] f_, t_, s_, st_, tlist_, slist_, options_ ] :>
  920.     mylaplace [ Cos[b] Cos[w t] f - Sin[b] Sin[w t] f, t, s,
  921.             st, tlist, slist, options ],
  922. mylaplace[ Cos[w_. t_] f_. CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  923.     Block [    {trans},
  924.         trans = mylaplace[f CStep[t], t, s, st, tlist, slist, options];
  925.         scaleL[addL[substituteForL[trans, s, s + I w],
  926.                 substituteForL[trans, s, s - I w] ],
  927.                1/2] ] /;
  928.     FreeQ[w, t],
  929.  
  930. (*          2.  multiplication by a sine function [Bracewell, 224]    *)
  931. mylaplace[ Sin[b_ + w_. t_] f_, t_, s_, st_, tlist_, slist_, options_ ] :>
  932.     mylaplace [ Sin[b] Cos[w t] f + Cos[b] Sin[w t] f, t,
  933.             s, st, tlist, slist, options ],
  934. mylaplace[ Sin[w_. t_] f_. CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  935.     Block [    {trans},
  936.         trans = mylaplace[f CStep[t], t, s, st, tlist, slist, options];
  937.         scaleL[subL[substituteForL[trans, s, s + I w],
  938.                 substituteForL[trans, s, s - I w] ],
  939.                I/2] ] /;
  940.     FreeQ[w, t],
  941.  
  942.  
  943. (*      H.  Derivative                        *)
  944. mylaplace[ Derivative[m_][f_][t_],
  945.             t_, s_, st_, tlist_, slist_, options_ ] :>
  946.     Block [ {dialogue},
  947.         dialogue = SameQ[ Replace[Dialogue, options], All ];
  948.         Assuming[ m > 0, dialogue ];
  949.         If [ dialogue && ! IntegerQ[m],
  950.              Print[ "and that ", m, " is an integer" ] ];
  951.         scaleL[ mylaplace[f[t], t, s, st, tlist, slist, options],
  952.             s^m ] ] /;
  953.     Implies[NumberQ[m], IntegerQ[m] && m > 0],    (* [Bracewell, 224] *)
  954.  
  955. mylaplace[ Derivative[m_][f_][t_] CStep[t_],
  956.             t_, s_, st_, tlist_, slist_, options_ ] :>
  957.     Block [ {dialogue},
  958.         dialogue = SameQ[ Replace[Dialogue, options], All ];
  959.         Assuming[ m > 0, dialogue ];
  960.         If [ dialogue && ! IntegerQ[m],
  961.              Print[ "and that ", m, " is an integer" ] ];
  962.         timeDerivative[mylaplace[f[t], t, s, st, tlist, slist, options],
  963.                    s, m, f[t], t] ] /;
  964.     Implies[NumberQ[m], IntegerQ[m] && m > 0],    (* [Bracewell, 224] *)
  965.  
  966.  
  967. (*      I.  Integral with time as upper limit                *)
  968. mylaplace[ Integrate[f_, terms1___, {a_, 0, t_}, terms2___],
  969.                t_, s_, st_, tlist_, slist_, options_ ] :>
  970.     Block [    {leftover, trans},
  971.         leftover = ToCollection[ terms1, terms2 ];
  972.         trans = mylaplace[f /. a -> t, t, s, st, tlist, slist, options];
  973.         If [ SameQ[{leftover}, {}],
  974.              scaleL [ trans, 1/s ],
  975.              scaleL [ trans, 1/s Integrate[f, leftover] ] ] ],
  976.                             (* [Churchill, 323] *)
  977.  
  978. mylaplace[ Integrate[f_, terms1___, t_, terms2___],
  979.        t_, s_, st_, tlist_, slist_, options_ ] :>
  980.     Block [    {leftover, trans},
  981.         leftover = ToCollection[ terms1, terms2 ];
  982.         trans = mylaplace[f, t, s, st, tlist, slist, options];
  983.         If [ SameQ[{leftover}, {}],
  984.              scaleL [ trans, 1/s ],
  985.              scaleL [ trans, 1/s Integrate[f, leftover] ] ] ],
  986.                             (* [Churchill, 323] *)
  987.  
  988.  
  989. (*      J.  Divide by powers of t                    *)
  990. mylaplace[ f_. t_^m_, t_, s_, st_, tlist_, slist_, options_ ] :>
  991.     integrateL[mylaplace[f t^(m+1), t, s, st, tlist, slist, options], s] /;
  992.     IntegerQ[m] && m < 0,                (*  [Nilsson, 516]  *)
  993.  
  994.  
  995. (*      K.  Conjugation:  L{ Conj f[-t] } --> Conj F(s)        *)
  996. mylaplace[ Conjugate[x_], t_, s_, st_, tlist_, slist_, options_ ] :>
  997.     conjL[mylaplace[x /. t -> -t, t, s, st, tlist, slist, options], s],
  998.                           (*  [Bracewell, 122]  *)
  999.  
  1000.  
  1001. (*      L.  Convert forms with Abs[a t] into two-sided transforms    *)
  1002. mylaplace[ f_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1003.     mylaplace[ absDialogue[f, t, options],
  1004.            t, s, st, tlist, slist, options ] /;
  1005.     ! FreeQ[ f, Abs[c_. t] ],
  1006.  
  1007.  
  1008.  
  1009.  
  1010. (*  V.    S T R U C T U R E S                        *)
  1011.  
  1012.  
  1013. (*      A.  Convolution                        *)
  1014. mylaplace[ CConvolve[t_][x1_, x2_, rest__], t_, s_,
  1015.         st_, tlist_, slist_, op_ ] :>
  1016.     multL[mylaplace[x1, t, s, st, tlist, slist, op],
  1017.           mylaplace[CConvolve[t][x2, rest], t, s,
  1018.             st, tlist, slist, op]],
  1019.  
  1020. mylaplace[ CConvolve[t_][x1_, x2_], t_, s_, st_, tlist_, slist_, op_ ] :>
  1021.     multL[mylaplace[x1, t, s, st, tlist, slist, op],
  1022.           mylaplace[x2, t, s, st, tlist, slist, op]],
  1023.  
  1024. mylaplace[ CConvolve[tconv_][x1_, x2_], t_, s_, st_, tlist_, slist_, op_ ] :>
  1025.     convolveL[ CConvolve[Complement[tconv, {t}]],
  1026.            mylaplace[x1, t, s, st, tlist, slist, op],
  1027.            mylaplace[x2, t, s, st, tlist, slist, op]] /;
  1028.     ListQ[tconv] && MemberQ[tconv, t],
  1029.  
  1030.  
  1031. (*      B.  Continuous FIR filters.                    *)
  1032. mylaplace[ CFIR[t_, h_List, Roots -> r_], t_, s_,
  1033.        st_, tlist_, slist_, options_ ] :>
  1034.     Block [ {i, numzeroes},
  1035.         If [ Length[h] != Length[r] + 1,
  1036.              Message[ CFIR::invalid ] ];
  1037.         numzeroes = Length[r];
  1038.         Transform[ Product[(s - r[[i]]), {i, 1, numzeroes}],
  1039.                -Infinity, Infinity ] ] /;
  1040.     Length[Dimensions[h]] == 1,
  1041.  
  1042. mylaplace[ CFIR[t_, h_, rest___], t_, s_, st_, tlist_, slist_, options_ ] :>
  1043.     Transform[ h . Table[s^k, {k, 0, Length[h] - 1}], -Infinity, Infinity ],
  1044.  
  1045. mylaplace[ CFIR[all__] [ x__ ], t_, s_, st_, tlist_, slist_, options_ ] :>
  1046.     multL[mylaplace[CFIR[all], t, s, st, tlist, slist, options],
  1047.           mylaplace[x, t, s, st, tlist, slist, options]],
  1048.  
  1049.  
  1050. (*      C.  One-dimensional analog IIR filters.            *)
  1051. (*          Region of convergence is incorrect            *)
  1052. mylaplace[ CIIR[t_, a_List, Roots -> r_], t_, s_,
  1053.         st_, tlist_, slist_, options_ ] :>
  1054.     Block [ {i, numpoles},
  1055.         If [ Length[a] != Length[r] + 1,
  1056.              Message[ CIIR::invalid ] ];
  1057.         numpoles = Length[r];
  1058.         rminus = Max[Re[r]];
  1059.         If [ ListQ[rminus], rminus = N[rminus] ];
  1060.         Transform[ 1 / Product[(s - r[[i]]), {i, 1, numpoles}],
  1061.                rminus, Infinity ] ] /;
  1062.     Length[Dimensions[a]] == 1,
  1063.  
  1064. mylaplace[ CIIR[t_, a_List], t_, s_, st_, tlist_, slist_, options_ ] :>
  1065.     Block [    {alist, denom, index, rminus, len},
  1066.         alist = ToList[a];
  1067.         len = Length[alist];
  1068.         denom = a[[1]] +
  1069.             Sum[ alist[[index]] s^(index - 1), {index, 2, len} ];
  1070.         rminus = Max[ Re[ GetRootList[denom, s] ] ];
  1071.         Transform[ 1 / denom, rminus, Infinity ] ] /;
  1072.     Length[Dimensions[a]] == 1,
  1073.  
  1074. mylaplace[ CIIR[all__] [x__], t_, s_, st_, tlist_, slist_, options_ ] :>
  1075.     multL[mylaplace[CIIR[all], t, s, st, tlist, slist, options],
  1076.           mylaplace[x, t, s, st, tlist, slist, options]],
  1077.  
  1078.  
  1079. (*      D.  Imaginary part of a sequence                *)
  1080. mylaplace[ Im[x_], t_, s_, st_, tlist_, slist_, options_ ] :>
  1081.     Block [    {trans},
  1082.         trans = mylaplace[x, t, s, st, tlist, slist, options];
  1083.         scaleL[ subL [ trans,
  1084.                    conjL[substituteForL[trans, s, -s], s] ],
  1085.             1/2 ] ],
  1086.  
  1087.  
  1088. (*      E.  Periodic sequence with period k                *)
  1089. mylaplace[ Periodic[a_, t_][f_] CStep[t_], t_, s_,
  1090.        st_, tlist_, slist_, options_ ] :>
  1091.     scaleL[mylaplace[f CPulse[a, t], t, s, st, tlist, slist, options],
  1092.            1 / ( 1 - Exp[- a s] ) ] /;
  1093.     FreeQ[a, t],
  1094.  
  1095.  
  1096. (*      F.  Real part of a function                    *)
  1097. mylaplace[ Re[x_], t_, s_, st_, tlist_, slist_, options_ ] :>
  1098.     Block [    {trans},
  1099.         trans = mylaplace[x, t, s, st, tlist, slist, options];
  1100.         scaleL [ addL [ trans,
  1101.                 conjL[substituteForL[trans, s, -s], s] ],
  1102.              1/2] ],
  1103.  
  1104.  
  1105. (*      G.  Reverse operator                        *)
  1106. mylaplace[ Rev[t_][x_], t_, s_, st_, tlist_, slist_, options_ ] :>
  1107.     antiCausalL[ mylaplace[ x, t, s, st, tlist, slist, options ], s ],
  1108.  
  1109.  
  1110. (*      H.  Shift operator                        *)
  1111. mylaplace[ Shift[m_, t_][x_], t_, s_, st_, tlist_, slist_, options_ ] :>
  1112.     scaleL[mylaplace[x, t, s, st, tlist, slist, options], Exp[- m s]] /;
  1113.     FreeQ[m, t],
  1114.  
  1115.  
  1116. (*      I.  Summation operator                    *)
  1117. mylaplace[ a_. Summation[i_, ib_, ie_, inc_][x_], t_, s_,
  1118.         st_, tlist_, slist_, options_ ] :>
  1119.     summationL[ mylaplace[a x, t, s, st, tlist, slist, options],
  1120.             Summation[i, ib, ie, inc] ],
  1121.  
  1122.  
  1123.  
  1124. (*  VI.   S T R A T E G I E S                        *)
  1125.  
  1126.  
  1127. (*      A.  Similarity --  f(at) <--> F(s/a) / |a|  [Bracewell, 224]    *)
  1128. (*         only perform this if a is neither 1 nor 0        *)
  1129. (*        must be done after operator rules because the        *)
  1130. (*        substitution t -> t/scale messes up operators        *)
  1131. mylaplace[ f_ CStep[t_], t_, s_, st_, tlist_, slist_, options_ ] :>
  1132.     Block [    {newf, scale},
  1133.         scale = ScalingFactor[f, t];
  1134.         newf = ( f /. t -> t/scale ) CStep[t];
  1135.         similarityL[ mylaplace[newf, t, s, st, tlist, slist, options],
  1136.                  s,
  1137.                  scale ] ] /;
  1138.     similarityQ[f, t],
  1139.  
  1140. mylaplace[ f_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1141.     Block [    {newf, scale},
  1142.         scale = ScalingFactor[f, t];
  1143.         newf = ( f /. t -> t/scale );
  1144.         similarityL[ mylaplace[newf, t, s, st, tlist, slist, options],
  1145.                  s,
  1146.                  scale ] ] /;
  1147.     similarityQ[f, t],
  1148.  
  1149.  
  1150. (*      B.  Place sub-expressions over a common denominator        *)
  1151. mylaplace[ f_. a_^(v_. t_ + k_.) / b_^(u_. t_ + l_.), t_, s_,
  1152.        st_, tlist_, slist_, options_ ] :>
  1153.     mylaplace[ (a^v/b^u)^t f a^k/b^l, t, s, st, tlist, slist, options ],
  1154.  
  1155.  
  1156. (*      C.  Expand out products like (t + 1) (t + 2) ...        *)
  1157. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1158.     mylaplace[ Distribute[x], t, s, st, tlist, slist, options ]  /;
  1159.     SameQ[Head[x], Times] && ! SameQ[Distribute[x], x],
  1160.  
  1161.  
  1162. (*      D.  Expands out all numerator in the expression        *)
  1163. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1164.     mylaplace [ Expand[x], t, s,
  1165.             SetStateField[st, expandfield, False],
  1166.             tlist, slist, options ]  /;
  1167.     GetStateField[st, expandfield] && ! SameQ[x, Expand[x]],
  1168.  
  1169.  
  1170. (*      E.  Expands out all numerators and denominators in expression    *)
  1171. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1172.     mylaplace[ ExpandAll[x], t, s,
  1173.            SetStateField[st, expandallfield, False],
  1174.            tlist, slist, options ] /;
  1175.     GetStateField[st, expandallfield] && ! SameQ[x, ExpandAll[x]],
  1176.  
  1177.  
  1178. (*      F.  Collect all terms in expression                *)
  1179. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1180.     mylaplace[ MyCollectAll[x, t], t, s,
  1181.            SetStateField[st, collectallfield, False],
  1182.            tlist, slist, options ] /;
  1183.     GetStateField[st, collectallfield],
  1184.  
  1185.  
  1186. (*      G.  Reduces all signal processing operations to their        *)
  1187. (*        mathematical forms                    *)
  1188. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1189.     Block [    {newx},
  1190.         newx = ToContinuous[ TheFunction[x] ];
  1191.         state = If [ SameQ[newx, x], st, initLstate[] ];
  1192.         state = SetStateField[state, thefunfield, False];
  1193.  
  1194.         mylaplace[ newx, t, s, state, tlist, slist, options ] ] /;
  1195.     GetStateField[st, thefunfield],
  1196.  
  1197.  
  1198. (*      H.  Apply definition if flag is set                *)
  1199. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1200.     Block [    {state, trans},
  1201.         state = SetStateField[st, definitionfield, False];
  1202.         trans = Integrate[x Exp[- s t], {t, -Infinity, Infinity}];
  1203.         If [ SameQ[Head[trans], Integrate],
  1204.              mylaplace[x, t, s, state, tlist, slist, options],
  1205.              definitionDialogue[ x CStep[t],
  1206.                      Transform[ trans, 0, Infinity ],
  1207.                      Integrate, options ] ] ] /;
  1208.     GetStateField[st, definitionfield] && Replace[Definition, options],
  1209.  
  1210.  
  1211. (*      I.  Two-sided transform if all else fails            *)
  1212. mylaplace[ x_, t_, s_, st_, tlist_, slist_, options_ ] :>
  1213.     Block [    {trans},
  1214.         trans = MyLaPlace[ x CStep[t] + x CStep[-t], t, s,
  1215.                    nullLstate[], tlist, slist, options ];
  1216.         If [ TrueQ[LForm[trans] && TheFunction[trans] != 0],
  1217.              MyMessage[Transform::twosided, trans, x],
  1218.              mylaplace[x, t, s, state, tlist, slist, options] ] ] /;
  1219.     GetStateField[st, stepfield]
  1220.  
  1221. }
  1222.  
  1223.  
  1224. (*  E N D     P A C K A G E   *)
  1225.  
  1226. End[]
  1227. EndPackage[]
  1228.  
  1229. If [ TrueQ[ $VersionNumber >= 2.0 ],
  1230.      On[ General::spell ];
  1231.      On[ General::spell1 ] ];
  1232.  
  1233.  
  1234. (*  A L I A S E S  *)
  1235.  
  1236. Unprotect[ LaPlaceTransform ]
  1237. LaPlaceTransform = LaPlace
  1238. LaPlaceTransform::usage = LaPlace::usage
  1239. Protect[ LaPlaceTransform ]
  1240.  
  1241.  
  1242. (*  H E L P     I N F O R M A T I O N  *)
  1243.  
  1244. Combine[ SPfunctions, { LaPlace } ]
  1245. Protect[ LaPlace ]
  1246.  
  1247.  
  1248. (*  E N D     M E S S A G E  *)
  1249.  
  1250. Print["The forward Laplace transform rule base LaPlace has been loaded."]
  1251. Null
  1252.